Boundary susceptibility in the open XX Z-chain 



Michael Bortzj and Jesko Sirkerf 

t Bergische Universitat Wuppertal, Theoretische Physik, 42097 Wuppertal, 
Germany 

| School of Physics, The University of New South Wales, Sydney 2052, Australia 
and 

Department of Physics and Astronomy, University of British Columbia, 
Vancouver, B.C., Canada V6T 1Z1 

Abstract. In the first part we calculate the boundary susceptibility \B m the 
open AXZ-chain at zero temperature and arbitrary magnetic field h by Bcthc 
ansatz. We present analytical results for the leading terms when \h\ -C a, where 
a is a known scale, and a numerical solution for the entire range of fields. In 
the second part we calculate susceptibility profiles near the boundary at finite 
temperature T numerically by using the density-matrix renormalization group for 
transfer matrices and analytically for T >C 1 by field theoretical methods. Finally 
we compare \B a t finite temperature with a low-temperature asymptotics which 
we obtain by combining our Bethe ansatz result with recent predictions from 
bosonization. 
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1. Introduction 

Even a single impurity can have a drastic effect on the low-energy properties of a 
one-dimensional interacting electron system. One of the simplest examples is an 
antiferromagnetic spin-1/2 chain with a non-magnetic impurity which cuts the chain 
and leads to a system with essentially free boundaries. Because translational invariance 
is broken the one-point correlation function (S z (r)) is no longer independent from the 
site index r and the local susceptibility %(r) acquires a nonzero alternating part [1]. 
Furthermore, the asymptotic behavior of correlation functions near such a boundary 
is no longer governed by the bulk critical exponents but instead by so called boundary 
or surface critical exponents [2, 3]. It is interesting to consider the case when the spin 
chain is not cut but instead one of the links is only slightly weaker. 



Here J > 0, J[ 2 > and we have allowed for an XXZ-type anisotropy which is 
described by the parameter 



N-l 



H = J 




(1) 



A =: COS7 with 0<A<1 (0<7< tt/2) . 



(2) 
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Using the Jordan- Wigner transformation one can also think about this system as a 
lattice model of spinless Fermions with a repulsive density-density interaction 



Therefore J[ = J — 5J with < 5 J <C J corresponds to a weakening of the hopping 
amplitude whereas J' 2 = J — 5 J gives a weakened density-density interaction along 
one bond. For all A both perturbations have the same scaling dimension x = K/2 
where K = tt/(-k — 7) is the Luttinger parameter [4]. Weakening the hopping or the 
interaction along one bond is therefore always relevantf. Assuming that the open 
chain presents the only stable fixed point one therefore expects that the physics at 
energies below Tk/J ~ {5J/J) X / X is governed by the open XXZ-chain [5]. This is 
the motivation to consider in the following only the open boundary condition (OBC) 



In an open XXZ-chain the free boundaries induce corrections of order 1/N to the 
bulk limit §. In particular a l/TV-term in the susceptibility is expected which we denote 
hereafter as boundary susceptibility XB{h, T). From the scaling arguments given before 
it follows that a long chain with a finite concentration of impurities is effectively cut 
into pieces of finite length. Measurements of the susceptibility on such a system will 
therefore reveal large contributions from the boundaries. This has inspired a lot of 
theoretical work to actually calculate these boundary corrections [7, 8, 9, 10, 11, 12]. 
Very recently the leading contributions to the boundary susceptibility for ft< 1 and 
Tel have been calculated by field theoretical methods [11, 12]. On the other hand it 
is also known that the XXZ-chain with OBC is exactly solvable by Bethe ansatz (BA) 
[13, 14]. For zero temperature, however, only the leading functional dependence on h 
for the isotropic case XB(h,T = 0) ~ l/h(\nh) 2 has been calculated so far [9, 10, 15]. 
For finite temperatures de Sa and Tsvelik [7] have applied the thermodynamic Bethe 
Ansatz (TBA) in the anisotropic case. Evaluating their TBA equations and comparing 
with a numerical solution (see section 4) we have found that their results are wrong 
for all anisotropics. Even the free Fermion case (see section 2) is not reproduced 
correctly and there is also disagreement with the field theoretical results by Fujimoto 
and Eggcrt [11] and Furusaki and Hikihara [12] at low temperatures. Frahm and 
Zvyagin [15] have treated the isotropic case with the same TBA technique. Although 
at least the functional form for low temperatures is correct, their results are not reliable 
for high temperatures [16]. This raises the question if the TBA is applicable at all for 
OBC or at least which modifications have to be incorporated compared to the well 
known case of periodic boundary conditions (PBC). 

Our paper is organized as follows: We start with the simple but instructive free 
Fermion case and establish results for the boundary susceptibility both as a function 
of T and h in section 2. In section 3 wc report the Bethe ansatz solution for T = and 
anisotropy < A < 1. We present analytical results for the boundary susceptibility at 
\h\ <C 1 and a numerical solution of the BA formulas for arbitrary h. In section 4 wc 
calculate susceptibility profiles near the boundary numerically by the density-matrix 

X For the free Fermion case K = 2 the perturbation becomes marginal. 

§ Note that in the periodic case no 1/N corrections exists and the first correction to the bulk limit 
is of order 1/N 2 and determines the central charge [6]. 




J[=J' 2 = 0. 
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renormalization group for transfer matrices (TMRG) and analytically by field theory 
methods. We also compare our numerical results for xs(^ = 0, T) with an analytical 
formula for T< 1 which we obtain by combining our BA results from section 3 with 
recent results from bosonization [11, 12]. The last section presents a summary and 
conclusions. 



2. Free spinless Fermions 



Here we want to consider the special case A = where cq. (3) describes noninteracting 
spinless Fermions. After Fourier transform the Hamiltonian takes the form 



N 

H = (J cos k n + h) (c\^ Ck n + h.c^] with k n = 



N + 1 



(4) 



where we have included a magnetic field h. Note that the only difference to PBC are 
the momenta k n which in this case would be given by k n — 2n/N. The susceptibility 
is easily obtained as 



^ (J cosfc„ + h) 



(5) 



and using the Euler-MacLaurin formula then yields 

1 



x(h,T) 
Xbulk 0, T) 
XB(h,T) 



Xbulk(h,T) + —XB(h,T) + O 



N' 



i r 
i 



i 



with 



cosh 



4ttT 
1 



cosh 



1 

2T 
1 

2T 



(J cos A: + h) 
(J cos A: + h) 



dk 
dk 



(6) 



8T 



cosh 



2T 



(J + h) 



1 

8T 



cosh 



2T 



(J-h) 



Therefore bulk and boundary susceptibility are identical at T — and given by 

1 1 



Xbulk(h,T = 0) = X B(h,T = Q) 



(7) 



Jtt y/i - (h/jy ' 

At finite temperatures Xbulk an< i Xb are different, however, the additional factors in 
Xb vanish exponentially for T — > so that they still share the same low-temperature 
asymptotics 



Xbulk(^-0,T^0) = XB (h = 0,T 



o) = ^ + ^ + o(n 



.(8) 



Fig. 1 shows the boundary and bulk susceptibilities for finite temperature at ft, = 
and as a function of h at T = (inset). In the next section we will discuss the 
BA solution for T = 0. We will see that a finite interaction between the Fermions, 
A / 0, has dramatic effects and Xbulk C 1 )^ 1 = 0) an d XB(h,T = 0) are no longer 
identical. For l/2<A<lwe find that XB{h 1 T = 0) even diverges for h — > whereas 
Xbulk {h = 0, T = 0) remains always finite. 




Figure 1. Bulk and boundary susceptibilities for free Fermions. Note that 
x(h,T = 0) diverges for h — > h c = J ■ 



3. The Bethe ansatz solution 



In this section we calculate ground-state properties of the model (1), i.e. T = 0. 
The Hamiltonian (1) with J[ = J' 2 = has been diagonalizcd both by coordinate 
and algebraic Bethe Ansatz [13, 14]. In the following, we refer to the algebraic Bethe 
Ansatz [14]. The eigenvalues E are parameterized by a set of M-many quantum 
numbers {Ai, . . . , Am}, 



E 



M 

E 



sin 2 7 



N — 1 



cos 7 



cos 7 



hS z 



(9) 



cosh(2Aj) 

S z = N/2-M. (10) 

Here, S z is the total z-component of the spin, h a magnetic field along the z-direction 
and we always assume in the following h > without loss of generality. The anisotropic 
exchange constant A is given by eq. (2) and the A& are determined by the following 
set of coupled algebraic equations 

0(A fc + i7/2) a(2A fc , 7 )a(A fc , tt/2 - 7 /2)a(A fc , tt/2 - 7 /2) 
<XA fe - i 7 /2) a(2A fc , -7)o(A fc , -tt/2 + 7/2)o(A fc) -tt/2 + 7 /2) 
(7M(A fc + ry)q M {-^k - 17) 



9M(Afe - i7)gM(-Afe + i7) ' 
with the definitions 



(11) 



M 



0(A) := sinh 2JV (A) , a(A, /z) := sinh(A + i/i) , 9m(A) := J^J sinh(A — Aj) . 

We first deal with the anisotropic case < 7 < 7r/2 and obtain equations for the 
susceptibility, which are solved analytically in the limit of small magnetic field. The 
isotropic case is treated afterwards. At the end of this section, we present numerical 
results for the susceptibility at arbitrary magnetic field. 
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3.1. Anisotropic case 

The solutions to (11) are periodic in the complex plane with period 27ri, so that we 
can focus on a strip parallel to the real axis with width 2ti\. Using arguments of 
analyticity, one sees that there are 2N + 3 + 2M roots in such a strip. So besides the 
M-many A& which yield the energy eigenvalues (9), there are 2N + 3 + M additional 
roots. Consider now the strip with Im Xk e] — n, tt] Vfc. We denote the roots in this 
set by {Ai, . . . , Am, X^\ . . . , )^ +M , 0, — i7r/2, i7r/2}. It is straightforward to verify 
that 0,±i7r/2 are roots. However, the algebraic Bethe Ansatz fails for these roots 
so that these solutions must be excluded. The roots are distributed symmetrically, 
both with respect to the real and imaginary axis. In this work we focus on the 
calculation of the ground state energy where M — N/2 and Xj e IR >0 Vj. Then there 
are N/2 roots A^ = — Xj, j = 1, . . . , N/2. A numerical evaluation of (11) shows that 

the remaining roots ^2n/2+i,...3N/2 ( A j=3iv/2+i,...5iv/2) nave imaginary part -\^/2 
(yy/2). The eigenvalues E are symmetrical in the Xj. We thus want to deal with 
the set {vi, . . . , v^} := {A^j 2 , . . . , x[ h \ Ai, . . . , Ajv/2}, whose elements are distributed 
symmetrically on the real axis w.r.t. the origin. From (11), we find that the Vj are 
the N real solutions to the equations 



ajvm, k/2 - 7/2) a(v m , -f/2) 
a(v m ,~ir/2 + 7/2) a(v m , -7/2) 



qN(Vm +17) 

qN{v m -h)' 



<P(v m + 17/ 2 ) 

4>{v m - 17/ 2 ) 

The remaining 2N solutions v^ 1 are identified as 2N = ^j=iv/2+i 5N/2- ^ n 

(12), the terms in brackets [• • •] are due to the OBC. These terms would be absent in 
the case of PBC. 

Our aim is to calculate the l/A^-contribution to the ground state energy per lattice 
site in the thermodynamical limit (TL). Like in the PBC-case [17, 18], we introduce 
the density of roots on the real axis, 

a v m +v m+1 V m - 1 +V m 

A v m ■■= 2 2 ' m = 2,...,N -1 

P + (v m ) := (13) 

where A Vm is the distance between two points on the left and on the right of the root 
v m , such that the left (right) point is situated midway between (v m ,v m+ i). 
From numerical studies the boundary values of p + are inferred, p+{v\) = p+(vn) = 0, 
with p + (vi,n)A Vi n = 1. Together with (13) it follows that 

N 1 

2 P+(v m )A Vm = - . (14) 

m— 1 

In the TL A^^ — > and p + (v m ) becomes a smooth function. Let us define the interval 
with non- vanishing density by [-B, B], i.e. v\ — » —B, — > B in the TL. Then sums 
over functions f(t>k) are transformed into integrals by 



E4/W= / mP+(x)dx- -^/(0) + 0(1/N% 
fc=l 



where the contribution /(0) is subtracted because the algebraic Bethe Ansatz fails 
at the origin. There are no further 0(l/iV)-terms because p + vanishes outside the 
integration boundaries by definition, 

p+ (v) = p+ (v)6(-v + B)6(B + v), 
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where 8(v) is the Heaviside-function. In order to find the continuum version of (12), 
it is convenient to define 

P(v) ■= P+{v)+P-{v) 
p-{v) := p(v)(6(v - B) + 6{-B - v)) . 
By taking the logarithmic derivative, the continuum version of (12) is obtained 

#(X, 7) + ^ 7) + #(X, 7T - 7) + #{X, 2 7 )] 

= p(x)+ [ d(x-y,2 7 )p + (y)dy, (15) 

J-B 

where 

„ . , 2isin7 d , sinh(a; + vy/2) . . 

27riz9(x,7 ) := — ' = — In— -7 r 1 ^. (16) 

cosh2x — cos 7 ax sinh(i — 17/2) 

Equation (15) is a linear integral equation with two unknowns, B and p. In a first 
step, (15) is solved for B = 00; in a second step, p+(x) is obtained depending on the 
parameter B and the dependence of B on the magnetic field h is calculated. We will 
see that B — 00 corresponds to h = 0, and a finite magnetic field h > induces a 
finite B < 00. Finally, the susceptibility x(/i) is deduced. This procedure was first 
used by Takahashi for PBC and is reviewed in [18]. 

Note that in deriving (15) the range of definition of the involved functions has been 
enlarged. All functions in (15) are defined on [—00, 00]; to calculate physical quantities 
(like the ground-state energy), however, we only need to know p + , defined on [-B, B}. 
Actually, it will be shown later that instead of dealing with p + , all quantities we are 
interested in can be expressed more conveniently by g+(x) := 9(x)p(v + B). The 
calculation of these functions is done by Fourier transformation, 

p(x) = —J_J(k)e- ikx dk. 

Let us first consider the case B = 00, where p = p + . It is straightforward to solve 
(15) in Fourier space, where 

~ _ sinh( 7 r/2-7/2)fc 
^ fc ' 7j ^ sinhTrfc/2 • 

We denote the solution of (15) for B = 00 by p n and find 

~m_ T m , 1 cosh 7 fc/4 cosh(7r/4 - 7 /2)fc 

po(fc) - s(k) + — cosh7fc/2 cosh( ^ _ i)k/A , (17) 



with 



s (^) : ~ o r — TTo' s ( x ) 



2cosh7fc/2' 27C0sh7rx/7 

Note that f™^ p(x)dx = 1/2 + 1/(2A), in agreement with (14). 

We now consider the case B < 00, i.e., a finite magnetic field. Let us derive the 
equation that determines g + . Using (17) we can rewrite (15) as 



p(x)=Po(x)+ n(x-y)p{y)dy (18) 

J\y\>B 

K ( X ) - J_ [°° sinh(^/2 - 7 )fc ikx 

[ ' ■ 2tt 2cosh 7 fc/2 sinh(7r - 7 )fc/2 6 [ ™> 
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We now introduce the functions 

p(x + B) =: g{x) = g+{x) + g_(x) 



g + (x) = 9(x)g(x), 
Then g(x) satisfies the equation 



.9-0) = 0(-x)g(x). 



(20) 



9{x) 



P CO POO 

= P a(x + B)+ / k(x - y)g{y) dy + / k{x + y + 2B)g{y) dy. (21) 
Jo Jo 

We seek a solution in the limit B 3> 0, which corresponds to \h\ <C a, where a 

is some finite scale which is calculated later. The driving term p (x + B) can be 

expanded in powers of exp [B]. Since (21) is linear in g and po, we make the ansatz 

9 = 9^ + 9^ + • ■ where superscripts denote increasing powers of exp [B]. Then 

gW(x) - [p (x + B)] (1) + k(x - y)gV>{y) dy 

Jo 

POO poo 

g^(x) = [p (x + B)} (n) + / K{x + y + 2B) g ^- 1 \y)dy+ / k{x - y)g {n \y) dy. 

Jo Jo 

Thus in each order, a linear integral equation of Wiencr-Hopf-type has to be solved. 

This technique is explained for example in [19, 20]. It relies on the factorization of 

the kernel k, 



1-k = 1/(G+G_), 



(22) 



where G + (G_) is analytical in the upper (lower) half plane and has asymptotics 
limi^i^oo G±(k) = 1. The functions G± are calculated in Appendix A. From (20), 
note that g + (g- ) is analytical in the upper (lower) half of the complex plane. Then 

g { +\k) = G + (k) [p (k)G-(k)e- ikB ] { ! ] (23) 



gf{k) = G + (k){[p«{k)G-{k)^ kB ] ( l ] 



where 



f±(k) :=± 



+ \^(k)g^(^k) G-(k)c- 2[kB G-(k) 
/(?) 



2vr 



dq 



(24) 



(25) 



, k — q ± ie 

is analytical in the upper (subscript +) or lower (subscript — ) half of the complex 
plane such that / = /+ + /_ . We will see later that it is sufficient to know g + . In 
this section we restrict ourselves to the calculation of g^ . The calculation of g^ is 
sketched in Appendix B. 



[. . in (23) is evaluated using (25), where only the pole nearest 
to the real axis is taken into account. Then we find 



The bracket 



9 ( +\k) = 



G+(k) 
1 



a 



+ 



2N 



k + in/j 



k + in/ 7 



+ 



k + i2tt/ (it — 7) 



-2tv B/(n-"f) 



for 7 7^ 7r/3 and 
g^(k) = G + (k) 



a 



-3B 



1 



k + i3ir 



2Nk + 3i 



Cl : Be~ 3B 



(26a) 



(266) 
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for 7 = 7r/3 where the constants are given by 

a = -G-{-mh) (27a) 

7 

\/2i ^ , . , s sin7r 2 /(47) 

°i = —G-{-mh) 7 9//A /,n 27& 

7 cos(7r 2 /(47) - 7r/4) v ' 

6i = — — tan7T7/(7r - 7) G_(-i27r/(7r - 7)) (27c) 

7T — 7 

1 8 

ci =i- 5 G_(-3i). (27d) 

For 7 = 7r/3, terms 0(£? exp [— £?]) occur in the boundary contribution which are 
absent for 7 = n/3. These terms are leading compared to those O (exp [— _£?]) so that 
these latter have been neglected for the boundary contribution in (266). 

We are now ready to compute s z := S z /N and e := E/N from (9), (10): 

s z = l/2-[ p+(x)dx + l/(2N) (28) 

J-B 

, „ J sin 7 f B n , . , . , J ( 2 — cos7\ 
e = -hs z ^± j ■&{x,- 1 )p+{x)dx + - (cos7+ jj^J- ( 29 ) 

We insert (18) into (28) to obtain 

= (0), (30) 

7T — 7 

which is an exact statement, including all orders . It is convenient to calculate 
e — e , where e := e(/i = 0). We use again (18) which yields 



, . 4j7rsin 7 f°° g + (x) 
e = -hs z H / ; v+v ^ dx 



7 J cosh(x + B)tt/^ 
hir (~(i) rn s ~(2)/ n x\ 87rJsin7 



7T — 7 



/ 

(#>(0) + g?>(0)) + [(^(iTT/7) +? (2) (W7)) 

- 5i 1) (3i7r/ 7 )e- 3 ' rB /^ + O (e" 3 ^/^ 2 ))] , (31) 

where in the last equation we restrict ourselves to the given orders. Now B is treated 
as a variational parameter and is determined in such a way that 

^(e-eo) = 0. (32) 

In this section we consider only the leading order in (31). Inserting (26 a), (30), (31) 
in (32), B is obtained as a function of h, 

B= -l]n- (33) 
7r a 

2nJ sin 7 G + (i7r/7) 

^ (7r ^ 7) ^7(or- (34) 

Thus a sets the scale for h. The restriction to the leading orders in cxp[— B] is 
equivalent to the leading orders in h in the limit \h\ <C a. 

One now makes use of (34) to determine s z (h) from (30), and therefrom x(/i) = 
ds z /dh. Inserting the explicit expressions for G± from Appendix A we find 
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The boundary contribution is given by 

f 7 sin7r 2 /(4 7 ) 




J(tt - 7)tt\/2 sin 7 cos(tt 2 /(47) - tt/4) 
7T7 1 r(7r/(7r — 7)) 



(Va) _(7r_37)/(7r_7) 



(356) 





(35c) 



and by 



1 /, 4/i \ 

" (,l) ^7^( ln ^ + 1 ) 



(35 d) 



for 7 = 7r/3. Note that the first term in (356), which is independent of magnetic field 
h, is the leading contribution for 7 > n/3 (pole closest to the real axis in (23)). For 
7 < 7r/3 the second term dominates and, in addition, this term is the next-leading 
contribution for 7 > tt/3 (second pole in (23)). For tt/7 < 7 < tt/3 the constant term 
represents the next-leading contribution, however, for 7 < ir/7 a pole at 67ri/(7r — 7) in 
(23) becomes second nearest to the real axis and gives the next-leading contribution 
(see Appendix B). 

The second term in (35 6) is in perfect agreement with the result obtained 
by conformal field theory and bosonization in [11, 12]. However, the first, field 
independent term has not been obtained before. Also the result (35rf) for the special 
case 7 = 7r/3 is new. Our results are in qualitative agreement with the TBA-work [7] 
at T — 0, where also a finite value of xb(T — 0, h = 0) for 7 > tt/3 (i.e. A < 1/2) and 
a divergent contribution with the same exponent as in (356) for 7 < tt/3 (A > 1/2) 
have been found. However, the coefficients in (356) differ from those in [7], due to 
an incorrect treatment of the Bcthe root at spectral parameter x — in [7] (c.f. the 
discussion in section 1). 

3.2. Isotropic case 

The isotropic case 7 = (i.e. A = 1) is treated in the same manner as the anisotropic 
case 7^0. For the bulk susceptibility, (35a), the limit 7^0 can be performed 
directly, yielding Xbulk = 1/Jtt 2 . For the boundary contribution this limit is more 
complicated and we describe the procedure in the following. 

First, we rescale (9) by Xj — ► -f\j. This is equivalent to substituting k — > fc/7 in 
Fourier space. Then 



Whereas the analyticity properties of the bulk contribution to (36) arc qualitatively 
the same as in (17), the boundary contribution shows, besides poles, a cut along the 
imaginary axis. In (23), the [. . .] + -brackct thus yields contributions 0(exp [—const. B]) 
from the poles, and algebraic contributions due to the cut. The exponential 
contributions are clearly sub-leading in comparison to the algebraic ones, so only the 



1 



2coshA:/2 




(36) 



Boundary susceptibility in the open XXZ-chain 



10 



latter are calculated in the following. Using eq. (A. 2) from Appendix A and explicit 
expressions for G±(k), we find (omitting the bulk contribution) 

' G + (0)(a 1 /B + a 2 (\nB)/B 2 + a 3 /B 2 ) 

+0((lnB)/B 3 ,l/5 3 ), k = 

\ ai G+{k)/{kB 2 ), k^O 

1 V2 i A i, ,„ , 

a.\ = , a 2 = ——^7 , «3 = — 7= — In 2 — - ln(2-7r) 

V2n 4tt 2 V2tt 2 V 2 y ' 

From (30), (31), (32) we obtain 

B =--ln- (37) 
7r a 

a -l G +(°) f o 8l 

° - 2kJG + ^) • (38) 

These equations are obtained by those from the anisotropic case, (33), (34), by scaling 
B — > 7_B and sending 7 — > afterwards. Carrying out the same steps which lead to 
(356), one finds the boundary contribution 



1 / 1 Ln|Ln/i/ft, 



4 W = - T TTTTT" + 777^777^ + >0 ( 39 ) 



ln/i//i 21n 2 /i/^o 

If 1 ln|ln/i//i | 1 \ / 1 \ 

XB(h) = -[ 5 + ' 1 01 ^ +o 5— 40) 

4V/iln 2 V^o h\n 3 h/h 2h\n 3 h/h J KhlifhJ 

h Q = a/V2 = Jny^Je . (41) 

The scale ho has been chosen such that in (39), no terms 0(ln~ 2 h) occur. The results 
(40), (41) agree with the TBA-work by Frahm et al. [15] for T = 0. Furthermore, 
agreement is found with [11, 12, 9], where scales which differ from ours (41) by a 
constant factor were used. 



3.3. Numerical evaluation 

To obtain results for the case when \h\ -^i a, (15) has to be solved numerically. 
For this purpose the bulk and boundary contributions to p{x) in equation (15) are 
treated separately. Both can be evaluated numerically for arbitrary values of B. 
The corresponding value for h is then derived from the minimum condition (32). 
This finally yields s z (h) and therefrom xW, c.f. eq. (28). The result for the bulk 
susceptibility is shown in fig. 2, together with the h = values (35 a). The boundary 
susceptibility is shown in fig. 3 (4) for A < 1/2 (A > 1/2). In both cases, the numerical 
solution confirms the analytical findings in the limit \h\ <^ a. 



4. Finite temperatures 

In the preceding section we calculated ground-state properties by making use of the 
intcgrability of the Hamiltonian (1) with J[ = J' 2 = 0. The next step would be to 
exploit integrability in order to calculate finite-temperature properties. However, the 
TBA seems to be problematic for systems with OBC as discussed in the introduction. 
The other available technique, namely the quantum-transfer-matrix-approach (QTM) 
[21], has not been applied to open systems so far. The problem to modify the TBA 
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H 0.3 




Figure 2. Bulk susceptibility from a numerical solution of eq. (15). The diamonds 
denote the h = values according to (35a). Note that the h = value is 
approached with infinite slope in the isotropic case due to logarithmic terms, 
cf. eq. (B.17). 




Figure 3. (a) Xs(^) f° r A = 0.0, ■ ■ ■ ,0.4. The stars denote the h = values 
according to eq. (356). (b) Comparison between the numerical solution (dots) of 
eq. (15) and the asymptotics (dashed lines) for \h\ <C a (eq. (356)). 




Figure 4. (a) XB(h) for A = 0.5, • ■ • , 1.0. (b) Comparison between the numerical 
solution (dots) of eq. (15) and the asymptotics (dashed lines) for \h\ <C a 
(cq. (356)). 



appropriately for OBC (if it is applicable at all) and the challenge to apply the QTM- 
method for OBC remain open issues for future research. Here, we use field theoretical 
arguments combined with a numerical study to discuss finite temperatures. 

First, we want to present a way different from section 3 to calculate the boundary 
susceptibility. Because translational invariance is broken in a system with OBC 
the one-point correlation function (S z (r)) is no longer a constant. The excess 
magnetization caused by the boundary can be defined as 

Mexc(r) = (S z (r)) OB C - ^PBC (42) 

where Mp-QQ is the magnetization per site in the system with PBC and r is the 
distance from the boundary. The local boundary susceptibility is then given by 
Xb{t) — dM exc (r)/dh and the total boundary susceptibility \b can be obtained 

by 

oo 

Xb = Xb{t) = xoBC - XPBC ■ (43) 

This means that we can calculate \b by considering only a local quantity which 
is particularly useful in numerical calculations where it is difficult to obtain the 
1/N contribution directly. Particularly suited for this purpose is the density- 
matrix renormalization group applied to transfer matrices (TMRG) because the 
thermodynamic limit is performed exactly. The idea of the TMRG is to express 
the partition function Z of a one-dimensional quantum model by that of an equivalent 
two-dimensional classical model which can be derived by the Trotter-Suzuki formula 
[22, 23]. For the classical model a suitable transfer matrix T can be defined which 
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allows for the calculation of all thermodynamic quantities in the thermodynamic limit 
by considering solely the largest eigenvalue of this transfer matrix. Details of the 
algorithm can be found in Refs. [24, 25, 26, 27]. The method has been extended to 
impurity problems in [28]. In particular, the local magnetization at a distance r from 
the boundary of a system with N sites is given by 



(S z (r)} 



N- 



£ n <*2|T"-iT|*£> 



(44) 



where 1^^.) ((^™|) are the right (left) cigcnstates of the transfer matrix T, T is 
a modified transfer matrix containing the broken bond and T(S Z ) is the transfer 
matrix with the operator S z included. Because the spectrum of T has a gap between 
the leading eigenvalue A and the next-leading eigenvalues, eq. (44) reduces in the 
thermodynamic limit to 



N 



lim (S z (r)) = 



(*2 / |T(5 z )T r - 1 T|*°j) 



(45) 

A5<*»|T|*») 

Therefore only the leading eigenvalue and the corresponding eigenvectors have to be 
known to calculate the local magnetization in the thermodynamic limit. Far away 
from the boundary (S z (r)} becomes a constant, the bulk magnetization 



m 



lim lim (S z {r)) 

r— >oo iV— »oo 



lim 



zj* L\ns z )T r - i \*i)(n\T\*°ii) 



(*°l\T(S')\*%) 



K(n\T\K) 



An 



(46) 



To obtain numerically the susceptibility profile xs( r ) at ft, = we calculate Mexc(?*) 
for small fields h/J^ 10~ 2 to 10~ 3 by using eqs. (45,46) and taking the numerical 
derivative. As an example we show in fig. 5 the susceptibility profile for A = 0.6 at 
various temperatures. 




Figure 5. Susceptibility profile for A = 0.6. The lines are a guide to the eye. 



At sufficiently low temperatures the susceptibility profile exhibits a maximum. 
This maximum gets shifted further away from the boundary as the temperature is 
lowered. The dependence of xs(r) on A at fixed temperature is shown in Fig. 6. Here 
the maximum becomes more pronounced with increasing A and is at the same time 
shifted further into the chain. 
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Figure 6. Susceptibility profiles for different A at T/J = 0.01387. The TMRG 
data are denoted by circles, the lines are the field theory results according to 
Eq. (51) and the dashed lines the field theory results shifted by some lattice 
spacings (see text below). 



Next we want to compare the numerics with field theory predictions. We start 
with the bulk two-point correlation function (S z (r)S z (0)) . The leading term in the 
long distance asymptotics of this function at zero temperature is known to be given 
by [29, 30, 31, 32, 33] 

(^(r)y(0))~A C ° s(2fc r f; + ^ (47) 

with an amplitude A and phase </>. The Fermi momentum is given by fcp = 7r(l±2m)/2 
and the scaling dimension by x — K/2. The usual mapping of the complex plane onto 
a semi-infinite cylinder then implies for small temperatures 

(g'(r)g'(0))~A ^ F r + 4>) K m 



( Jk, S i n h il> 



Using Cardy's relation between 2n-point functions in the bulk and n-point functions 
near a surface [2] one can now directly obtain the magnetization near the boundary 

(S'(r))~A ™( 2k * r + +l /2 . (49) 



Note that although the critical exponent is only half the exponent appearing in the 
two-point bulk correlation function both decay exponentially with exactly the same 
correlation length £ = v s /(2irxT). With the known result for the bulk susceptibility 
XbulkC 1 = 0) = K/2irv s we obtain 

<S'(r)> ^i H)rSm(flr/ ;;! ^ h« land 
(^sinh^> N 
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AK (-l) r r 
XB(r)\ h =o — -jot ■ (50) 



This is the leading contribution to the boundary susceptibility. Note that eq. (50) 
agrees for the special case A = 1 with the result given in [1]. The amplitude depends 
only on the operator product expansion of S z and is given by A — \J~AJ2 where A z 
has been derived by Lukyanov and Terras [33] (see eq. (4.3)). However, this alternating 
term does not contribute when calculating \b by integrating over all lattice sites. The 
leading non-oscillating contribution has already been obtained by Fujimoto and Eggert 
[11] and Furusaki and Hikihara [12]. Including this term we find for the susceptibility 
profile 



Xb(t)\ 



h=0 



a z k {-iy 



(^sinh^>) 



K/2 



4K 2 \ r 2 

+ -T2-7 72K (51) 

Vs (^^h^>) 

where the amplitude A is given in [33, 12]. This field theory result is shown as straight 
lines in Fig. 6 in comparison to the numerics. For all A the shape of the curves 
agrees well with the numerical results. However, especially at larger A the height of 
the maximum is overestimated and there is also a shift by a few lattice sites. When 
we shift the theoretical curves by an appropriate amount of lattice spacings (dashed 
lines in Fig. 6) we see that the predicted exponential decay for larger distances agrees 
perfectly with the numerical data. 

First, we should note that we cannot expect that the field theoretical treatment 
yields reliable results for short distances. In addition, the next-leading alternating 
terms in the asymptotic expansion of the bulk two-point correlation function will 
become more and more important as A — * 1 [33]. This makes our approximation 
to take only the leading term (47) into account apparently worse for larger A. In 
fact shifting the field theoretical result is equivalent to taking contributions with 
larger scaling dimensions into account. Our observation that the shift increases with 
increasing A is therefore consistent with the increasing importance of next-leading 
terms. We also want to mention that a similar shift has been observed before for the 
isotropic case [1]. 

Finally we want to discuss the total boundary susceptibility \b at finite 
temperature. To calculate it numerically we have in principle to sum xb(t) over 
all lattice sites. However, at the lowest considered temperature the correlation length 
£ < 50 and it is sufficient to take the sum over the first 200 sites around the boundary. 
The results for A = 0.1 — 0.4 are shown in Fig. 7 and for A = 0.6 — 0.9 in Fig. 8. For 
the low-temperature asymptotics we already know from our Bethe ansatz calculations 
that there is a temperature and magnetic field independent term which dominates for 
A < 0.5. In addition there is a temperature dependent contribution which stems from 
the integral over all lattice sites of the non-oscillating term in eq. (51) and dominates 
for A > 0.5. This integral has already been calculated in [11, 12], yielding the leading 
temperature-dependent contribution to \b{T). However, the constant term cannot 
be calculated within the field-theoretical framework. Our knowledge of this term from 
the BA calculations in the previous section allows us to obtain a low-temperature 
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0-1 T/J 



0-1 T/J 



Figure 7. Comparison between TMRG (circles) and the low-tempcraturc 
asymptotics (lines) according to (52) for A = 0.1, • ■ ■ , 0.4. 




Figure 8. Same as fig. 7 with A = 0.6, • • ■ , 0.9. 
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asymptotics which is valid for all < A < 1. This asymptotics is given by 

K-l sinf,^- 
XB (T,h = 0)= ' ' 



JV2tt sin(K - 1) cos f 

N CT(^)r(3-2^)(^-2^(X)) [2t I T\ 2K - 3 
+ Av^2-1/K)T(2-K) \ v ) l& > 

where ^>'(x) — d 2 lnr(a;)/da; 2 . The first term is our BA result, eq. (356), where 
we have substituted 7 = ir(K — 1)/K. The second term is the leading temperature 
dependent contribution taken from [11, 12]. As in the case h 7^ 0, T = 0, we expect 
that these are the leading and the sub-leading contributions for A < 0.9. For larger A 
we expect that other temperature dependent terms will become more important than 
the constant term and we have to omit it in this case to stay consistent. Note that 
the leading temperature dependent term has the same exponent 2K — 3 as the leading 
magnetic field dependent term in eq. (356). This is consistent with scaling arguments. 
The lines in figs. 7,8 denote the low-T asymptotics according to this formula and are in 
excellent agreement with the numerics except for A = 0.9 where the asymptotic limit 
at the lowest considered T is not yet reached. However, this is expected due to the 
afore mentioned other temperature dependent terms which will become sub-leading 
for A > 0.9 and even equally important as the term with exponent 2K — 3 for A — > 1, 
yielding finally a logarithmic dependence on temperature [11, 12]. 



5. Conclusions 



In the first part we have calculated the boundary contribution to the magnetic 
susceptibility of the XXZ-chain with OBC at zero temperature and finite magnetic 
field by BA. For small magnetic fields and 7 < tt/3 (A > 1/2) the BA result for the 
leading divergent term agrees with the field theoretical analysis [11, 12]. For 7 > tt/3 
(A < 1/2) a field independent term is dominating. This term has not been obtained 
before. We also derived for the first time the leading term for the special case 7 = tt/3 
(A = 1/2). In addition we have presented a numerical solution of the BA equations 
for arbitrary field h. We used the numerics for a verification of our analytical results 
for \h\ < a. 

In the second part we have calculated numerically susceptibility profiles near the 
boundary by the TMRG method and compared these results with a low-temperature 
asymptotics which we obtained by field theoretical methods. Apart from a shift by 
a few lattice sites we have found good agreement. By combining a temperature 
and magnetic field independent term, which we obtained by BA and which is 
dominating for A < 1/2, with the leading temperature dependent term, which has 
been calculated in [11, 12] and dominates for A > 1/2, we have obtained a low- 
temperature asymptotics for xb(T) which is valid for all < A < 1. Numerically, 
Xb(T) has been obtained by a summation of xb(t, T) over a sufficient number of sites 
around the boundary. At low-temperatures excellent agreement with the analytical 
formula was found. The remaining challenge is to calculate the finite-temperature 
properties analytically by using the integrability of the model, either by TBA or by 
the QTM-method. 

We acknowledge very helpful discussions with Andreas Klumpcr, Ian Affleck and 
Frank Gohmann. 
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-iak 



7 

— ln(-7r/7 — 1) — ln(l — 7/71") 

7T 



Appendix A. Factorization of the kernel 

Let us first carry out the factorization (22) for the anisotropic case 
1 sinh7rfc/2 

G+(k) G-(k) = 2cosh 7 A:/2 sinh(7r - 7 )fc/2 

G-{k) =G+(-k). 
Using properties of the T function, we find 

r (k)- fit T r(l-ifc/2) 

+ ( )-V \* - 7) r(1/2 _ i7jfc /( 27r )) r(l - ijfc(7r - 7)/(2tt)) 
1 

° =2 

where a is determined such that liniifci^oo G±(k) = 1. 

As already mentioned in section 3.2, the isotropic limit is realized by scaling 
k — > fe/ 7 , thus for 7 = 0, 

G+(fc)G_(fc) 2coshfc/2■ 1 ■ ' 

By making use of 

-T = !'""=- 1 1 "*"'*). < A ' 2 > 

the exponential in (A.l) is factorized in functions analytical in the upper and lower 
half planes with 

, (-ifcl-^A 271 ") . , 

^W-^ rwa I ■*/(*)) • 

1 ln(27r) 
° 2^ 2?r 

Appendix B. Next-leading orders 

Our focus here is on the anisotropic case; we comment on the isotropic case at the 
end of this appendix. 

The calculation of the next-leading order, i.e. of in (24), is technically more 
involved than the leading order g^ x \ because there are two contributions in (24). The 
calculation is done following the same steps as in section 3, so that we merely give the 
results here. 

The [• • -] + -brackets in (24) are evaluated using the integral representation (25). 
Now, the pole next-nearest to the real axis is taken into account in the first summand 
in (24). The second term already contains a factor exp [— 2\kB\, so that the pole next 
to the real axis yields the leading contribution. Thus we find for 7 ^ tt/3 

g { +\k) = G + {k) 



\ V k + i3r 



1 

~2N 



+ 



+ "0,2 \ c -37rB/ 7 _j Q 0,3 c -( 7r / 7 +47r/(7r- 7 ))B 

i37r/ 7 k + m/jj k + \2ty/{tt — 7) 

Ql » 1 + ai - 2 \ C-37TB/7 + a h3 e _( 7r/7+4w /( 7 r- 7 ))B 

k + i37r/7 k + m/jj k + i27r/(7r — 7) 

61,1 



k + i6ir/ (n — 7) 
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ao,i 



k + i2tt/ (tt — 7) 

7 



., — 677/ (77 — 7)73 



n,2 



k + iir/j 



-2(7r/7+7r/(7T-7))B 



7 



1 n ^3 

ao,2 = - — tan — G_ 

Z77T 2-f 



. ft 

7 



1 *. 7r 7 ^ 
a o 3 = —7 r tan G_ 

7r(yT + 7) 7T — 7 

. 2 sin7r/4 sin37r 2 /(47) 
7 cos(37r 2 /(47) + 7r/4) 

°i . 71-2 ^2 

°1,2 = T^tan— G_ 
2/T 27 



ai,i = 



-i- ) G^ 

7 



. 2tt 

L 

7T — 7 



7 



ai,3 = 
h,i = i 



017 



ti(tt + 7) 
2 



tan 



7T7 
7T — 7 



G 2 _ -i 



7T — 7 



tan i2!LG_ 

7T — 7 
^2 



27T 

7T — 7 
6/T 



-1- 



_ bl(7T-7) TT- 2 

61 2 - —7 — - — r tan — G_ 

7T(/T + 7) Z7 

fel3 = ^ tan ^2_G 2 

47T 7T — 7 

and for 7 = 7r/3 
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7T — 7 
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. 2tt 
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7T — 7 
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G_(-9i), ci, 



3 ci 



2 Gi(-3i). 

7T° 7T^ 7T Z 

This expression for <?2 is inserted into (31), where we now have to keep all the indicated 
terms. || Then B as a function of h is derived. In section 3, we found that this 
relationship is the same both for the boundary and for the bulk in the leading order. 
This is no longer true when next-leading terms are considered. For the bulk we obtain 



47/(71—7) 



a 0:2 G-(-i37r/ 7 ) 

A\ = h . , A 2 

«o G_ (-171-/7) 

and for the boundary 
ae -^/7 = h(l + A l 

l+2 7 /(7r- 7 ) 



^2 



tt - 7 «0,3 

27 a 

^\ 47/(^-7) 
a 



+ 



1 ao,4 
3 a 



- I ln- 

a 



(B.12) 



- 7— 7r/3 y 



+ B 1 



+B 2 - 



Cl,2 



/i ci G_(-9i) 



— l+67/(7T — f) 



- ) ln- 

a 



(B.13) 



3ai \a J a 3 G_(— 3i) 

_ 7=77/3 f 

|| It is indeed sufficient to restrict the expansion of the 1/ cosh(x + _B)-factor in (31) to the first two 
orders. The next term would involve the exponent 577/7. Comparing with the largest exponent in 
(B.l), 5-77/7 > 77/7 + 477/(77 — 7) for 7 < 77/2. However, 7 = 77/2 is allowed, since in this case, all 
coefficients except ao,i, ai,i vanish. 
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A 2 = 



it - 7 ai j3 



ffli,2 G_(-i37r/7) 

ai G-{— iir/j) ' 27 ai 

B = 2(tt - 7) 6i£ ^ = 7T + 7 61,2 G_(-i37r/7) 61 

7T + 7 01 ' 2 2(77 — 7) °i i7r/7) ai 



In (B.12), (B.13) and in the following, for 7 = 7r/3 the only nextleading terms are 
those labeled by [. . .] 7=7r y 3 . Combining these equations with (30), one finds 



s bulkW - 



tt(it — 7) 



G-(-Wt) 
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- tan — G a 
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(B.15) 



7— tt/3 ^ 

From these expressions, \ can be obtained. Let us consider the special case of free 
Fermions, i.e., 7 = it/2. Eq. (19) implies n = so that G+ = G_ = 1. Furthermore, 
from (27a)-(27d), (B.l)-(B.IO), the only non-vanishing coefficients are 



a = 1- , ai = 1- , a .i 

TT 7T 

so that 



-1- , ai,i 

7T 



-2B 
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-i-j 

7T 
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A; + 6i 



(l + i) + o(.-). 
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Note that this is the expansion of p (x + B) in Fourier space. Equations (B.12), (B.13) 
are equivalent in the free-Fermion case and yield 



c- 2B = -|l 
a 



Finally, the sum of (B.14), (B.15) can be simplified to 



xW = 



(B.16) 



where we have set a = 2 J for 7 = 7r/2. Equation (B.16) is in agreement with the 
exact result (7) within the first two orders. 

As far as the isotropic case 7 = is concerned, note that (39), (40) include already 
next-leading terms for the boundary magnetization and susceptibility. Logarithmic 
corrections to the finite bulk susceptibility eq. (35a) with 7 = 0, have been calculated 
by BA-techniques for PBC [34]: 



XbulkO) = 



Jtt 2 



1 + 



1 



lnln(/i //i) 



2\n{h /h) 4(ln(/i //i)) 2 



+ o(ln- 2 (/.))(B.17) 



with ho — ae 1 / 8 7r 1 / 4 . The scale h has been determined such that no terms 
O (ln~ 2 h) appear in (B.17). 
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